Equilibrium phase behavior of polydisperse hard spheres 
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We calculate the phase behavior of hard spheres with size polydispersity, using accurate free energy 
expressions for the fluid and solid phases. Cloud and shadow curves, which determine the onset of 
phase coexistence, are found exactly by the moment free energy method, but we also compute the 
complete phase diagram, taking full account of fractionation effects. In contrast to earlier, simplified 
treatments we find no point of equal concentration between fluid and solid or re-entrant melting 
at higher densities. Rather, the fluid cloud curve continues to the largest polydispersity that we 
study (14%); from the equilibrium phase behavior a terminal polydispersity can thus only be defined 
for the solid, where we find it to be around 7%. At sufficiently large polydispersity, fractionation 
into several solid phases can occur, consistent with previous approximate calculations; we find in 
addition that coexistence of several solids with a fluid phase is also possible. 

PACS numbers: 82.70.Dd, 64.10.+h, 82.70.-y, 05.20.-y 



During the past few decades, a great deal of effort has 
been devoted to studies of the phase behavior of spheri- 
cal particles, and in particular of the freezing transition, 
where the particles arrange themselves into a crystal with 
long-range translational order. The simplest system for 
studying this transition is one where the particles act as 
hard spheres, exhibiting no interaction except for an in- 
finite repulsion on overlap. This scenario can be realized 
experimentally, using e.g., colloidal latex particles steri- 
cally stabilized by a polymer coating Q. Hard spheres 
constitute a purely entropic system; the internal energy 
U vanishes, and F = —TS. Phase transitions are thus en- 
tropically driven; nevertheless, monodisperse (i.e., iden- 
tically sized) hard spheres exhibit a freezing transition, 
where a fluid with a volume fraction of 'd> « 50% coexists 
with a crystalline solid with cf> ss 55% 0. 

For colloidal hard spheres, there is inevitably a spread 
in the particle diameters a, which are effectively continu- 
ously distributed within some interval. The width of the 
diameter distribution can be characterized by a polydis- 
persity parameter S, defined as the standard deviation of 
the size distribution divided by its mean. 

The effect of polydispersity on the phase behavior of 
hard spheres has been investigated by experiments Q, , 
computer simulations @, 0, U@ , density functional the- 
ories 0, HI, and simplified analytical theories 0, 0, ITfl 



mmOlEi; Ref. 13 has a more detailed bibliography 
of earlier work. These studies have revealed that, com- 
pared to the monodisperse case, polydispersity causes 
several qualitatively new phenomena. First, it is intu- 
itively clear [jj that significant diameter polydispersity 
should destabilize the crystal phase, because it is dif- 
ficult to accommodate a range of diameters in a lat- 
tice structure. Experiments indeed show that crystal- 
lization is suppressed above a terminal polydispersity of 
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St w 12% 0,0. Theoretical work suggested that this 
arises from a progressive narrowing of the fluid-solid co- 
existence region with increasing 5, with the phase bound- 
aries meeting at 8t 0, llOf in a point of equal concentra- 
tion Bartlett and Warren |13| also found re-entrant 
melting on the high-density side of this point: for 6 just 
below St, they predicted that compressing a crystal could 
transform it back into a fluid, as sketched in the inset of 
Fig. H below. However, none of these theoretical studies 
fully accounted for fractionation 15], i.e., the fact that 
coexisting phases generally have different diameter dis- 
tributions; in fluid-solid coexistence, one typically finds 
that the solid contains a higher proportion of the larger 
particles. Beyond the resulting difference in mean diam- 
eter, fractionation implies that coexisting phases can also 
have different polydispersities S. Indeed, numerical simu- 
lations that allow for fractionation show that a solid with 
a narrow size distribution can coexist with an essentially 
arbitrarily polydisperse fluid 0, , suggesting that the 
concept of a terminal polydispersity is useful only for the 
solid but not for the fluid. Fractionation has also been 
predicted to lead to solid-solid coexistence |lllll2| . where 
a broad diameter distribution is split into a number of 
narrower solid fractions. This occurs because the loss of 
entropy of mixing is outweighed by the better packing, 
and therefore higher entropy, of crystals with narrow size 
distribution; accordingly, as the overall polydispersity of 
the system grows, the number of coexisting solids is pre- 
dicted to increase. 

Previous work as described above leaves open a num- 
ber of questions. The drastic and differing approxima- 
tions for size fractionation used in the studies of re- 
entrant melting and solid-solid coexistence 0, 0, 
leave the relative importance of these two phenomena 
unclear. In [13j fractionation was allowed, but coexisting 
phases were implicitly constrained to have the same 5; 
calculations that account fully for fractionation remain 
restricted to highly simplified van der Waals free ener- 
gies 0] ■ Numerical simulations have been carried out at 
constant chemical potential distribution 4. 16]; in con- 
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trast to the experimental situation, the overall particle 
size distribution can then change dramatically across the 
phase diagram, limiting the applicability of the results. 

Our goal in this letter is to calculate the equilibrium 
phase behavior of polydisperse hard spheres on the basis 
of accurate free energy expressions, taking full account of 
fractionation and going beyond previous work on fluid- 
solid and solid-solid coexistence. The experimentally ob- 
served behavior of hard sphere colloids will of course also 
depend on non-equilibrium effects, e.g., the presence of 
a kinetic glass transition [TtI ]. anomalously large nucle- 
ation barriers or the growth kinetics of polydisperse 
crystals ^tj. Nevertheless, the equilibrium phase be- 
havior needs to be understood as a baseline from which 
non-equilibrium effects can be properly attributed. Also, 
more of the equilibrium behavior may be observable un- 
der microgravity conditions, where the glass transition is 
shifted to higher densities or even absent pol |. 

Our calculations will show that the fluid cloud curve, 
which locates the onset of phase coexistence coming from 
low density, continues to large polydispersities 5: the 
point of equal concentration found in [13j disappears to- 
gether with the predicted re-entrant melting. Instead of 
returning to a single- phase fluid at high volume fractions, 
the system splits into two or more fractionated solids, 
consistent with the simplified calculations of jnj; coex- 
istence of several solids with a fluid phase appears as a 
new feature. 

In general, the total free energy (density) of a poly- 
disperse system consists of an ideal and an excess part 
(/ex)' In units where k^,T = 1, 

f = J dap(a)[]np(a)-l)+f eK (1) 

Here p(a) is the density distribution, i.e., p(a)da is the 
number density of particles with diameters between a 
and a + da. Equilibrium requires equality of the chemical 
potentials fi(a) = Sf/5p(a) and of the pressure II = — /+ 
f da p(a) p(o~) among all coexisting phases a — I... P. 
Particle conservation adds the condition that, if phase 
a occupies a fraction t>W of the system volume, then 
Y,a v(a) P ia) (v) = /O (0) (ct), where p (0) (a) is the overall or 
"parent" density distribution. 

For the fluid, the most accurate free energy approx- 
imation available at present is the BMCSL general- 
ization |2lj, |22j of the monodisperse Carnahan-Starling 
equation of state. This is truncatable in the sense that 
the excess free energy only depends on the four moments 
Pi = [da o % p{a) (i — ... 3) of the density distribu- 
tion [23j; po is the total number density, (ir/6)ps = 4> 
the volume fraction, and pi/po = a and P2/P0 = c 2 give 
the mean and mea n-sq uare diameter. For the crystalline 
solid, Bartlett 0, l24j assumed that the same truncat- 
able structure holds; an approximate excess free energy 
(depending only on the same pi) can then be derived 
from simulation results |25j for bidisperse hard spheres. 
Implicit in the use of data from 25] is the assumption 



that the crystal has a substitutionally disordered f.c.c. 
structure. 

We adopt the BMCSL and Bartlett free energies for 
our calculation; the appropriate branch for a given p(a) 
is selected by taking the minimum of the fluid and solid 
free energies. Since the excess free energies depend only 
on the pi, the excess chemical potentials p ex (a) take the 
form 

Mex(c) = Sf cx /Sp(a) = po + /iier + /i 2 cr 2 + pzo- 3 (2) 

For the solid, Bartlett [24| derived [1q and p% from 
the small and large a limits of the Widom insertion 
principle |26|. However, because of the approximate 
character of the excess free energy, p, cx (a) then does 
not obey the thermodynamic consistency requirement 
5 p ex (a) / S p(a' ) = Sp ex (a')/Sp(a). To avoid this, we as- 
sign all excess chemical potentials by explicitly carrying 
out the differentiation in 

Our computational a ppr oach is based on the moment 
free energy method |23, l28l |29| , which maps the full free 
energy Q, with its dependence on all details of p{a) 
through the ideal part, onto a moment free energy de- 
pending only on the moments pi. For truncatable free 
energies this locates exactly the cloud points, i.e., the 
onset of phase separation coming from either a single- 
phase fluid or solid, as well as the properties of the co- 
existing "shadow" phases that appear there. Inside the 
coexistence region, one in principle needs to solve a set of 
highly coupled nonlinear equations |l5j] and the moment 
free energy method gives only approximate results. How- 
ever, by retaining extra moments with adaptively chosen 
weight functions [2jJ |3(], , increasingly accurate solu- 
tions can be obtained by iteration. Using these as initial 
points, we are then able to find full solutions of the exact 
phase equilibrium equations. Care is taken to check that 
solutions are globally stable, i.e., that no phase split of 
lower free energy exists [2|j. We are able to calculate 
coexistence of up to P = 5 phases, which so far has been 
possible only for much simpler free energies depending 
on a single density moment (see e.g., |2{| ). 

Below we present results for a symmetric triangular 
parent density distribution, i.e., p^(cr) increasing lin- 
early from zero for a € [1 — w, 1] and decreasing linearly 
for a 6 [1,1 + w], with w — ^/66. The mean diame- 
ter of 1 fixes our length unit. Other distributions could 
be considered, but for the moderate values of S of inter- 
est here one expects them to give qualitatively similar 
results, based on the intuition that for narrow size dis- 
tributions 5 is the key parameter controlling the phase 
behavior ||J. 

Fig. ^ shows our results for the cloud and shadow 
curves. The fluid cloud curve continues throughout the 
whole range of polydispersities that we can investigate: 
even at S = 14%, a hard sphere fluid will eventually split 
off a solid on compression. Fractionation is key here; as 
indicated in Fig. ^ the coexisting shadow solid always 
has a smaller polydispersity, with 5 never rising above 
6%. This fractionation effect prevents the convergence of 
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FIG. 1: Cloud (thick) and shadow (thin) curves, plotted as 
polydispersity S versus volume fraction <f>; dashed lines link 
sample cloud-shadow pairs. The fluid (F) cloud curve con- 
tinues up to the largest 5 that we study. The solid (S) cloud 
curve has two branches, with onset of F-S and S-S coexistence 
at low and high volume fractions, respectively. Inset: Sketch 
of the phase diagram of fl3| , showing re-entrant melting and 
the point of equal concentration. 

the solid and fluid phase boundaries, along with the re- 
sulting re-entrant melting p| (Fig.nj-inset). These find- 
ings are in qualitative accord with numerical simulations 
for the simpler case of fixed chemical potentials |a, [l6| . 
In particular, the terminal polydispersity St cannot be 
defined as the point beyond which a fluid at equilibrium 
will no longer phase separate; St only makes sense as the 
maximum polydispersity at which a single solid phase 
can exist. As in [5j we also find that the coexisting fluid 
always has a lower volume fraction than the solid, along 
with (not shown) a lower mean diameter. 

Coming from the single-phase solid, decreasing den- 
sity at low polydispersities leads to conventional fluid- 
solid phase separation. At higher 5, however, the solid 
cloud curve acquires a second branch at higher densi- 
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FIG. 2: Full phase diagram for polydisperse hard spheres with 
a triangular size distribution. In each region the nature of 
the phase(s) coexisting at equilibrium is indicated (F: fluid, 
S: solid). Dashed line: best guess for the phase boundary in 
the region where our numerical data become unreliable. 



ties. This is broadly analogous to the re-entrant phase 
boundary found in |13| . but with the crucial difference 
that the system phase separates into two solids rather 
than a solid and a fluid. The two branches meet at a 
triple point. Here the solid cloud phase coexists with 
two shadow phases, one fluid and one solid, as marked 
by the squares in Fig.Q] From Fig.^the triple point, at 
St ~ 7%, also gives the terminal polydispersity beyond 
which solids with triangular diameter distribution are un- 
stable against phase separation. As explained, other dis- 
tributions should give similar values of St- 

In Fig. |21 we show the full phase diagram for our tri- 
angular parent distribution. In each region the nature of 
the phase(s) coexisting at equilibrium is indicated. The 
cloud curves of Fig.^reappear as the boundaries between 
single-phase regions and areas of fluid-solid or solid-solid 
coexistence. Starting from the latter and increasing den- 
sity or S, fractionation into multiple solids occurs. The 
overall shape of the phase boundaries in this region is 
in good qualitative agreement with the approximate cal- 
culations of ^lj. However, the coexisting solids do not 
necessarily split the diameter range evenly among them- 
selves as assumed in ^l] ; see the sample plot in Fig. 
of the normalized diameter distributions n(a) = p(a)/po 
of four coexisting solids. In fact, plotting S vs cf> for all 
coexisting solids across the phase diagram, we find points 
that cluster very closely around the high-density branch 
of the solid cloud curve in Fig. ^ Coexisting solids with 
lower volume fraction 4> thus tend to have higher poly- 
dispersity 6, as in the example in Fig. |3 this conclusion 
is intuitively appealing since higher compression should 
disfavor a polydisperse crystalline packing. 
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FIG. 3: Normalized diameter distribution of four coexisting 
solid phases obtained from a parent with (<j), S) — (63.0%, 8%). 
From left to right, the solids have volume fractions and 
polydispersities (60.1%, 5.4%), (62.9%, 4.6%), (64.6%, 4.0%), 
(66.3%, 3.6%). 

Note that in Fig. |21 at larger 5 than we can tackle 
numerically, coexistence of P > 4 solids would be ex- 
pected since each individual solid can only tolerate a 
finite amount of polydispersity. However, from Fig. [21 
such phase splits would occur at increasing densities and 
eventually be limited by the physical maximum volume 
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FIG. 4: Normalized diameter distributions for F-S-S-S phase 
coexistence obtained from a parent with (</>, 5) = (60.3%, 8%). 



fraction </> max ~ 74%. Also, at higher S more compli- 
cated single-phase crystal structures, with different lat- 
tice sites occupied preferentially by (say) smaller and 
larger spheres, could appear and compete with the sub- 
stitutionally disordered solids we consider. 

Finally, a new feature of the phase diagram in Fig. 
is the coexistence of a fluid with multiple solids. The 
triple point on the solid cloud curves already indicated 
the existence of a three-phase F-S-S region; as in the 
case of solid-solid phase splits, more solid phases then 
appear with increasing 5. Fig. 0] shows again that the 
fractionation behavior is non-trivial: while the coexisting 
fluid is enriched in the smaller particles as expected, it 
also contains "left over" large spheres that did not fit 
comfortably into the solid phases and thus ends up having 
a larger polydispersity (10.4%) than the parent (8%). 

In conclusion, we have calculated the phase behavior 



of polydisperse hard spheres, using accurate free ener- 
gies for the fluid and solid phases and solving exactly the 
resulting equilibrium conditions. Fluid-solid coexistence 
has been identified for fluids with polydispersities up to 
S = 14%. This shows clearly that the experimentally 
observed suppression of crystallization above 5 = 12% 
is a non-equilibrium effect, probably caused by increased 
nucleation barriers at large S |18|. For the solid, a termi- 
nal polydispersity remains well-defined as the maximal 
value beyond which instability to phase separation sets 
in; for triangular diameter distributions this turns out to 
be 8t pa 7%. Instead of the re-entrant melting predicted 
in an approximate treatment of fractionation effects |13| , 
we find that sufficiently polydisperse solids split into two 
fractionated solids on compression. At higher volume 
fractions and polydispersities, multiple solids can coex- 
ist; coexistence of a fluid with several solids appears as 
a new feature. Fractionation effects are nontrivial, with 
solids splitting the diameter range unevenly among them 
and coexisting fluids sometimes having larger polydisper- 
sities than the parent. 

Overall, our calculated phase diagram unites, clarifies 
and extends the previous separate predictions of polydis- 
persity effects on fluid-solid coexistence and solid-solid 
fractionation. Numerical simulations may offer the best 
avenue for testing our predictions but will need to be car- 
ried out at fixed parent size distribution [22j to detect the 
complex fractionation phenomena we find. For the future 
it would be exciting to unify our predictions with those 
for fluid-fluid demixing, but this will be very challenging 
since the latter only occurs at polydispersities 6 of order 
100% [13,03, far outside the range studied here. 

It is a pleasure to thank Paul Bartlett for providing his 
code for the solid free energy. Financial support through 
EPSRC grant GR/R52121/01 is acknowledged. 
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